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Summary. 

Material indentation studies, in which a probe is brought into controlled physical contact with 
an experimental sample, have long been a primary means by which scientists characterize 
the mechanical properties of materials. More recently the advent of atomic force microscopy 
which operates on the same fundamental principle, has in turn revolutionized the nanoscale 
analysis of soft biomaterials such as cells and tissues. This paper addresses the inferential 
problems associated with material indentation and atomic force microscopy through a frame- 
work for the changepoint analysis of pre- and post-contact data that is applicable to experi- 
ments across a variety of physical scales. A hierarchical Bayesian model is proposed to ac- 
count for experimentally observed changepoint smoothness constraints and measurement er- 
ror variability, with efficient Monte Carlo methods developed and employed to realize inference 
via posterior sampling for parameters such as Young's modulus, a key quantifier of material 
stiffness. These results are the first to provide the materials science community with rigor- 
ous inference procedures and uncertainty quantification, via optimized and fully automated 
high-throughput algorithms, implemented as the publicly available software package BayesCP. 
To demonstrate the consistent accuracy and wide applicability of this approach, results are 
shown for a variety of data sets from both macro- and micro-materials experiments — including 
silicone, neurons, and red blood cells — conducted by the authors and others. 

Keywords: Changepoint detection; Constrained switching regressions; Hierarchical Bayesian 
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1. Introduction 

This article develops a hierarchical Bayesian approach for contact point determination in 
material indentation studies and atomic force microscopy (AFM). Contemporary applica- 
tions in materials science and biomechanics range from analyzing the response of novel 
nanomaterials to deformation (Wong et al., 1997), to characterizing disease through me- 
chanical properties of cells, tissues, and organs (Costa, 2004). Experimental procedures and 
analyses, however, remain broadly similar across these different material types and physical 
scales (Gouldstone et al., 2007), with the scientific aim in all cases being to characterize 
how a given material sample deforms in response to the application of an external force. 

As illustrated in Figure 1 overleaf, indentation experiments employ a probe (or can- 
tilever arm, in the case of AFM) to apply a controlled force to the material sample. This 
indenting probe displaces the sample while concurrently measuring resistive force, with the 
resultant force-position data used to infer material properties such as stiffness (in analogy 

^Address for correspondence: Patrick J. Wolfe, Statistics and Information Sciences Laboratory, 
Harvard University, Oxford Street, Cambridge, MA 02138-2901, USA 
E-mail: wolfe@stat.harvard.edu; software URL: http://sisl.seas.harvard.edu/BayesCRhtml 



D. Rudoy, S. G. Yuen, R. D. Howe and R J. Wolfe 



IX 



Force Sensor 



Material Sample 




Fig. 1. Diagram of a macro-scale indentation experiment in which a spherical probe, attached to a 

force sensor, indents a material sample and deforms it by a distance 5. In this hypothetical example, 
a net change of 1.2 Newtons in resistive force is consequently observed. 



to compressing a spring in order to experimentally determine its spring constant). Prior 
to subsequent data analysis, a key technical problem is to determine precisely the moment 
at which the probe comes into contact with the material. Sample preparation techniques 
and sizes frequently preclude the direct measurement of this contact point, and hence its 
inference from indenter force-position data forms the subject of this article. 

At present, practitioners across fields lack an agreed-upon standard for contact point 
determination; a variety of ad hoc data pre-processing methods are used, including even 
simple visual inspection (Lin et al., 2007a). Nevertheless, it is well recognized that precise 
contact point determination is necessary to accurately infer material properties in AFM 
indentation experiments (Crick and Yin, 2007). For example, Dimitriadis et al. (2002) 
show that for small displacements of thin films, estimation errors on the order of 5 nm for a 
2.7 /Ltm sample can cause an increase of nearly 200% in the estimated Young's modulus — the 
principal quantifier of material stiffness. When soft materials such as cells are studied at 
microscopic scales, for example to determine biomechanical disease markers (Costa, 2004), 
the need for robust and repeatable AFM analyses becomes even greater (Lin et al., 2007a). 

In this article, we present the first formulation of the contact point determination task as 
a statistical changepoint problem, and subsequently employ a switching regressions model to 
infer Young's modulus. Section 2 summarizes the basic principles of material indentation, 
showing that the resultant force-displacement curves are often well described by low-order 
polynomials. Section 3 introduces a corresponding family of Bayesian models designed to 
address a wide range of experimental conditions, with specialized Markov chain Monte Carlo 
samplers for inference developed in Section 4. Following validation of the proposed inference 
procedures in Section 5, they are employed in Section 6 to infer material properties of mouse 
neurons and human red blood cells from AFM force-position data. The article concludes in 
Section 7 with a discussion of promising methodological and practical extensions. 



2. Material indentation 

2. 1. Indentation experiments and data 

Indentation experiments proceed by carefully moving a probe from an initial noncontact 
position into a material sample, as shown in Figure 1, while measuring the resistive force 
at some prescribed temporal sampling rate. After a small deformation is made, the probe 
retracts to its initial position; during this stage the resistive force decreases with every 
subsequent measurement. At the conclusion of each such experiment, two force-position 
curves are produced, examples of which are shown in Figure 2. In this article we consider 
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Distance from initial probe position (mm) Distance from initial probe position (jim) 

(a) Silicone indentation data (b) Red blood cell indentation data 



Fig. 2. Examples of force-position curves from which material stiffness properties are to be inferred, 

with force measurements during indentation and subsequent probe retraction shown in blacl< and 
grey, respectively, and subsets of the indentation data shown near the presumed contact points. 
Experiments were performed using a mechanical arm for the soft silicone sample, and an atomic 
force microscope for the red blood cell sample; note the differences in physical scale and noise level. 



only the forward- indentation curves, as is standard practice (Lin et al., 2007a), though the 

methods wc present arc extensible to retraction data whenever suitable models arc available. 

Despite significant differences in physical scale and noise level, the curves of Figure 2 
(and indeed, most indentation data sets) feature several common characteristics. In the 
pre-contact region, force response appears linear in the position of the indenter; drift due 
to experimental conditions is often present, yielding a non-zero slope as in Figure 2(a). The 
post-contact data appear well modeled by low-order polynomial functions of the correspond- 
ing displacement, and indeed, the key assumption of such experiments is that, conditioned 
upon knowledge of the geometry of the indenting probe, the relationship between the degree 
of material deformation and the measured resistive force depends in a known way on the 
material stiffness. This is analogous to the case in which an ideal spring is compressed a 
specified distance (5 by a force of known magnitude; a measure of the spring's stiffness is 
given by the spring constant k, which may be calculated using Hooke's law: F = —k6. 



2.2. Contact mechanics and the Hertzian model 

Indentation data such as those shown in Figure 2 are typically acquired for linear- elastic 
materials. This implies not only that the material instantaneously returns to its original 
shape following the cessation of an external force, but also that the relationship between the 
applied stress (force per unit area) and the resultant strain (deformation per unit length) is 
linear. This ratio is known as Young's modulus — the primary quantifier of material stiffness 
introduced above — and is reported in units of Pascals. 

In small-deformation experiments to infer Yoimg's modulus E for linear-elastic materials, 
both the indenter geometry and the measured resistive force come into play, by way of the 
so-called Hertzian model (Lin et al., 2007a). Specifically, the relationship between the 
sample deformation depth 6 and the measured resistive force F takes the form 



F(xE-6^, 



(1) 
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where the constant of proportionaUty and the fixed parameter /? depend on the geometry 
of the indenter tip in a known way. Examples include indentation by a sharp pyramid with 
tip angle 2(p, or a sphere of radius R, whereupon (1) takes the following forms: 

F=^^^^E.S' (pyramid) or F = ^^^^ E ■ 5'/' (sphere), (2) 

with V a known, dimensionless quantity termed Poisson's ratio. A subsequent fitting of 
the Hertzian model of (1) to experimental data thus allows one to obtain an estimate for 
Young's modulus E once the post-contact region has been identified. (Below we retain the 
standard practitioner notation {E, F, R, 6, v, (f>), as distinct from other variables to follow.) 

2. 3. Pre- and post-contact data regimes 

Observe that the Hertzian model of (1) posits a relationship between force and indentation 
depth, whereas the data of Figure 2 are seen to be a function of the position of the indenter. 
The Hertzian model thus describes the underlying physics of the post-contact stage of a 
typical indentation experiment, while the measured data also comprise a pre-contact stage. 
As we detail below, the union of these two regimes is well described by a switching regressions 
scenario, with force measurements prior to contact typically linear in the position of the 
probe. The corresponding intercept represents the equilibrium tare required to achieve zero 
net force, while the slope captures constant-velocity drift in force measurements that can 
arise in a variety of experimental settings. 

In both the pre- and post-contact data regimes, it is standard to assume independence of 
measurement errors — a reasonable assumption for all practically achievable force sampling 
rates. Errors prior to contact arise due to force sensor vibration in the experimental medium, 
thermal variations, and other effects — whereas after contact they also depend on interactions 
between the probe and the sample such as frictional forces. Consequently, while error 
variances in these regimes can be expected to differ in practice (analysis of AFM data from 
a red blood cell, shown in Figure 6 and described later, reveals one such example), their 
relative magnitudes are not known a priori. As a final consideration, uncertainty in the 
reported position of the indenting probe is typically several orders of magnitude smaller 
than the distance between consecutive force sampling points, and can safely be disregarded. 

As noted by Crick and Yin (2007), all of the sources of variability mentioned above 
can lead to large relative errors in resistive force measurements near the contact point; 
Figure 2(b) illustrates a typical scenario. This makes manual identification of the contact 
point difficult in many cases, and motivates the model-based approach that we now describe. 

3. Bayesian changepoint model 

Having outlined the basic principles of material indentation, we now proceed to formulate a 
Bayesian model for force-position data that encompasses both pre- and post-contact regimes. 
In light of the discussion above, the corresponding task of contact point determination is 
recognizable as a changepoint estimation problem in the context of switching regressions. 
In turn, the Hertzian model of (1) implies that an estimate of Young's modulus can be 
obtained as a linear function of the leading post-contact regression coefficient. Control 
over experimental conditions implies that each indentation data set contains precisely one 
changepoint, thus obviating any need to estimate the presence or number of contact points. 
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During an indentation experiment, the indentcr moves continuously through a sequence 
of n positions x = {xi,X2, - ■ ■ and records a force measurement yi at each position 
Xi, resulting in a sequence of force measurements y = (yi, j/2, • • • , As we shall later 
treat models for soft material indentation, in which the pre- and post-contact curves are 
constrained to be continuous at the regression changepoint, we begin by introducing a 
continuous parameter 7 e (l,n) denoting the contact point index, with the corresponding 
contact point at which the indenter first contacts the sample denoted by x^ e {xi,Xn)- 



3. 1. Data likelihood for indentation experiments 

We adopt a classical switching regressions scenario for our model, in which y is assumed 
to be a polynomial function of known degree di in position x prior to contact, and of 
known degree ^2 in deformation depth S = x — x^ after contact with the sample is made. 
This formulation encompasses the Hertzian model of (1) if fractional powers are allowed; 
however, for clarity of presentation we consider ^2 to be an integer unless otherwise noted. 
Letting p = di+d2 + 2 denote the number of regression coefficients in our model, and with 
n the number of data points, the corresponding design matrix is hence of dimension nx p. 
We subsequently employ the subscript 7 to denote any quantity that depends on 7, and 
index via subscripts 1 and 2 the pre- and post-contact regression regimes. 

We denote the regression coefficients by /3i e M'^i+-^ and /32 G R^'^-i-i^ with design 
matrices Xi^j,X2^-y defined as follows, for [7J the largest integer less than or equal to 7: 



/I 

1 
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X2 



\1 X 



hi 



X 
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1 
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The observed data y may likewise be partitioned into pre- and post-contact vectors 

Via = ?y2, • • • , V\_-i\ ) and ^2,7 = (j/L7j+i' ?>L7j+2' • • • ' ' 



(3) 



and, following our discussion in Section 2.3 regarding the noise characteristics typical of 
indentation experiments, we assume independent and Normally distributed additive errors, 
with unknown variances a\,a2- Thus for 1 < i < n, we have that is distributed as follows: 



Vi 



[Ar(Xi,^/3i,a2) ifl<i<L7j, 
\M{X2a^2.ol) if L7j+l<^<ri. 



(4) 



Note that the statistical model of (4) is consistent with the Hertzian mechanics model 
of (1), resulting in a post-contact force-response curve that is a power of the displacement 
b = X — x^ oi the material sample, rather than the position x of the indenter. However, 
when (^2 is an integer, the cocifficicnts of these two polynomials are related by a simple linear 
transformation. Consider, for instance, a quadratic curve in 6 given by f{6) = ao+ai6+a2S'^ . 
One may rewrite f{6) as a quadratic polynomial in x as follows: 



Oo + ai{x — x^) + a2{x — x^Y = 60 + bix + b2X^ 



(5) 



where 62 = 02, &i = ai — a2Xj, and bo = ao — aix^ 

X2a to be reformulated directly in terms of indenter position x, such that it no longer 



a2xi^. This transformation enables 
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depends continuously on x^, in contrast to (3). Transformations akin to (5) do not apply, 
however, when the d2 is a fraction, as in the case of (2) for a spherical indenter, or when a 
continuity constraint is enforced at the changepoint; we detail such cases below. 



3.2. A general parametric Bayesian model for material indentation 

The likelihood of (4), together with the presence of genuine prior information dictated by 
the underlying physics of material indentation experiments, suggests a natural hierarchical 
Bayesian model. In contrast to the semi-conjugate approach taken by Carlin et al. (1992), we 
detail below a fully conjugate model, as this allows for analytical simplifications that we have 
observed to be important in practice. Integrating out nuisance parameters improves not 
only the mixing of the chains underlying the resultant algorithms and inferential procedures, 
but also their computational tractability when data sizes grow large. 

We specify prior distributions for all model parameters, including the contact point 
index 7 G (l,n), the pre- and post-contact regression coefBcients /3i and /32, and the error 
variances and cr|. For i £ {1,2}, we then assume that f3i ^ A/'(/Xi, crf A~^), with a 
{di + 1) X {di + 1) diagonal positive definite matrix and /iti G R'''+^. A standard inverse- 
Gamma conj ugate prior IQ (ao , 60 ) is adopted for both variances ai and trl . Finally, recalling 
that a single regression changepoint is always assumed within the data record, wc employ 
a uniform prior distribution on the interval (l,n) for the contact point index 7. In certain 
experimental settings, whereupon the initial position of the indenter is known to be at least 
a certain distance from the sample, an informative prior distribution may well be available. 

Because of sensitivity to hyperparameters, we follow standard practice and adopt hyper- 
priors for increased model robustness. A Gamma prior is assumed on 60 so that bo ^ Q{k, rj), 
however, we determined from simulations that the posterior estimators considered were not 
sensitive to the prior parameters Hi and Aj, and so did not employ an additional level of 
hyperprior hierarchy for the regression coefficients. For notational convenience, we let Ifc 
denote the 1-vector of dimension k whose entries are all equal to one and the zero matrix 
of appropriate dimension, and define the variables y, (3, fi, X^, A, S^, and S as follows: 



= diagKlL^j,cT2l„_L^j) e R"^", S = diag(a?ld,+i,c7|ld,+i) e W^'p. 

The posterior probability distribution of the model parameters (7, /3, df , ctI, 60), conditioned 
on the observations y and the fixed model parameters ip = (/i, A, oq, k, rf), is then: 

p{'y,f3,cri,ff2,bo\y;tp) ocp(y|/3,o-?,(Ti,7)p(^|fTf,ffi;/x, A)p(cTi|6o;ao)p(o-i|6o;oo)p(6o;K,??)p(7) 

■ (|S^| |S|)-i exp [-1 {{y-X^fSy {y-X^O) + {(3- n)' E'^A (/3 - m)}] ■ (6) 

To confirm robustness, we also studied the effect of replacing the diagonal prior covari- 

ancc A for the pre- and post-contact regression coefficients (3i and (32 with an appropriately 
adapted 5-prior (Zellner, 1986) such that /3i|pi,7 ~ M{iJ,i,a^pf{Xl^^Xi^j)~^), with a 
scale parameter to which we ascribed a diffuse inverse-Gamma hyperprior. We observed 
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no measurable effect of this change in priors on the resulting inference further confirming 
the insensitivity of the adopted model to the prior covariance of the regression coefficients. 
Moreover, efficient sampling from the conditional distribution of pf is precluded by its de- 
pendence on the contact point index 7, reducing the overall efficacy of this approach in the 
Markov chain Monte Carlo approaches to posterior sampling described in Section 4 below. 

3.3. Smoothness constraints at the changepoint 

It is often the case that force-position curves are continuous at the contact point x^. Es- 
pecially for soft materials such as the red blood cells we consider in Section 6, it is to 
be expected that the change in the force measurement is smooth, and a continuity con- 
straint can serve to regularize the solution in cases where many different fits will have high 
likelihood. Imposing smoothness constraints dates back to at least Hudson (1966) who con- 
sidered this constraint in deriving maximum likelihood estimators for switching regressions. 
More recently, Stephens (1994) used it in a hierarchical Bayesian setting. 

In our setting, according to the likelihood of (4), a continuity constraint on the pre- and 
post-contact force-position curves at x = implies that 

/3io + fhixj + ■■■+ Pid, x^' = P20 , (7) 

where f3i = (/3io,/?ii, . • . ,/3idi)' and (32 = (/32o,/32i, • • • ,024^)' denote the vectors of pre- 
and post-contact regression coefficients, respectively. Higher-order smoothness can also 
be imposed: we say that the force-position curve is s times continuously differentiable at 
x-y if the sth-order derivatives of the pre- and post-contact curves meet at x^, with (7) 
corresponding to the case s = 0. On the other hand, if ^2,7 were a function of the position 
X, rather than the displacement x — x.y, then the continuity constraint would become: 

di d2 

J2Pu^\ = J^P2jXl,. (8) 
i=0 j=0 

Either continuity constraint implies that the likelihood function is nonlinear in the contact 
point x^; enforcing more degrees of smoothness at the changepoint serves to exacerbate 
the nonlinearity and makes the design of efficient inference methods increasingly difficult. 
Stephens (1994) imposed (8) in a Bayesian switching regressions setting and proposed a 
rejection sampling step within a Gibbs sampler to address the resultant nonlinearity. Later, 
in Section 4, we describe a more efficient approach that can be applied when either (7) 
or (8) (or higher-order analogues) are enforced. 

3.4. Changepoint estimation and contact point determination in the literature 

As demonstrated above, inference for material indentation data is well matched to classical 
statistical frameworks for changepoint estimation. Independent of the specifics of our con- 
tact point problem, the last half century has seen a vast body of work in this area. Sequential 
and fixed-sample-sizc varieties have been considered from both classical and Bayesian view- 
points, with numerous parametric and nonparametric approaches proposed. We refer the 
interested reader to several excellent surveys, including those by Hinkley et al. (1980), Za- 
cks (1983), Wolfe and Schechtman (1984), Carlin et al. (1992), and Lai (1995). Some of 
the earliest work on maximum likelihood estimation of a single changepoint between two 
polynomial regimes was done by Quandt (1958) and Robison (1964). 
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Historically, Clicrnoff and Zacks (1964) were among the first to consider a parametric 
Bayesian approach to changepoint estimation. Changepoints arising specifically in linear 
models have been treated by many authors, including Bacon and Watts (1971), Ferreira 
(1975), Smith (1975), Choy and Broemling (1980), Smith (1980), and Mcnzcfrickc (1981). 
The introduction of Markov chain Monte Carlo methods has led to more sophisticated hi- 
erarchical Bayesian models for changepoint problems, beginning with the semi- conjugate 
approach taken by Carlin ct al. (1992), in which the prior variance of the regression co- 
efi[icients is left unsealed by the noise variance. Advances in trans-dimensional simulation 
methods have rekindled interest in multiple changepoint problems, as discussed by Stephens 
(1994), Pimskaya et al. (2002) and Fcarnhcad (2006), among others. 

In the context of material indentation, however, existing approaches to contact point de- 
termination do not make use of changepoint estimation methodology. In fact, current meth- 
ods arc error-prone and labor intensive — even consisting of visual inspection and manual 
thresholding (Lin et al., 2007a). However, as described earlier, the many sources of variabil- 
ity in indentation data imply that one cannot always simply proceed "by eye." Moreover, 
in the context of atomic force microscopy, most experiments aiming to characterize cell 
stiffness, for example, employ multiple, repeated indentations at different spatial locations. 
These requirements have motivated a more recent desire for effective, high-throughput au- 
tomated techniques, as detailed in (Lin ct al., 2007a). 

Interpreted in a statistical context, the procedures thus far adopted by practitioners 
fall under the general category of likelihood fitting. Rotsch et al. (1999) suggest simply 
to take two points in the post-contact data and solve for E and 7 using the appropriate 
Hertzian model; however, the resultant estimate of Young's modulus depends strongly on 
the indentation depth of the selected points (Costa, 2004). Costa et al. (2006) propose to 
minimize the mean-squared error of a linear pre-contact and quadratic post-contact fit to 
the indentation data, though under the assumption of equal pre- and post-contact variances. 

None of the existing approaches adopted by practitioners, however, provides any means 
of quantifying uncertainty in changepoint estimation — an important consideration in prac- 
tice, since measurement errors can be large relative to the reaction force of the probed mate- 
rial, and consequently may result in poor point estimates (Crick and Yin, 2007). Moreover, 
such approaches fail to capture necessary physical constraints of the material indentation 
problem, such as the smoothness constraints described in Section 3.3 above. Such short- 
comings provide strong motivation for the hierarchical model developed above, as well as 
the robust and automated fitting procedures we describe next. 

4. Posterior inference via lUlarl^ov chain IVIonte Carlo 

The hierarchical Bayesian modeling framework introduced above features a large number of 
unknowns, with constraints on certain parameters precluding closed-form expressions for the 
marginal posteriors of interest. These considerations suggest a simulation-based approach 
to inference; indeed, it is by now standard to use Markov chain Monte Carlo methods to 
draw samples from the posterior in such cases. Though widely available software packages 
for Gibbs sampling arc adequate for inference in certain hierarchical Bayesian settings, the 
complexity of the conditional distributions we obtain here (after imposing constraints and 
integrating out parameters whenever possible) necessitates explicit algorithmic derivations 
on a casc-by-case basis. To this end, we build upon the approaches of Carlin ct al. (1992) 
and Stephens (1994), and employ Metropolis- within- Gibbs techniques to draw samples from 
the posterior of (6) as well as under the smoothness constraints of Section 3.3. 
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4. 1. Metropolized Gibbs Samplers and Variance Reduction 

The selection of conjugate priors in our model allows nuisance parameters to be integrated 
out, in order to reduce the variance of the resultant estimators. Following standard manip- 
ulations, we marginalize over the pre- and post-contact regression coefficients /3i and (32, 
respectively. This yields the following marginal posterior probability distribution: 

2(ao-l)^-6o/<T| ^-1 60/,, 



P('y,cri,(72,bo\y;tp) oca-^ 'e ^ ■ ^ " 'e ^ ■ bg 



e 



■ (|S^| |E| l^^l)-^ exp {-i {y'S-'y + m'S'V - b^'A-'b^)} , (9) 

where = X^'S-^X,^ + S-^A e is block-diagonal and = AS" V + ^t'S^^?/ G 

]gj)x 1 rpjjg marginal posterior of (9) factors into a Gamma density in 60, and inverse-Gamma 
densities in and erf by way of the following partitions of A^ and b^: 



The cixprc'Ksions of (9) and (10) lead to thc^ following Gibbs sampk^: 



Algorithm 1 Gibbs sampler for changepoint estimation 

(a) Draw 7 ~ p(7|(Ti , (t|, 60, y; 1/') according to (9); 

(b) Draw ^Xg{ao + ^[l},bo + |(yi,7'yi,7 +Ati'Mi - bi^j' A^l^bi^j)); 

(c) Draw trf ~ J^(ao |(n- L7j),^'o + 5(z/2,7'z/2,7 + - &2,7'^^i^b2,7)); 

(d) Draw /)„ - t/f/v. + a^^ + rr^^) ■ 



To simulate from the conditional distribution of 7, we employ as a Metropolis-within- 
Gibbs step a mixture of a local random walk move with an independent Metropolis step 
in which the proposal density is a pointwisc evaluation of (9) on the grid 1.2,...,n of 
indenter location indices. It is also possible to integrate out both noise variances (or the 
hyperparameter bo)- In this case, additional Metropolis steps are required, as the resulting 
conditional density of is nonstandard. Our simulation studies confirm that these variants 
exhibit less Monte Carlo variation than a Gibbs sampler based on the full posterior of (6). 



4.2. Posterior inference in the presence of smoothness constraints 

Bearing in mind the underlying physics of soft materials, a sufficiently high sampling rate 

of the force sensor relative to the speed of the indenter may yield data consistent with a 
smoothness assumption of a given order s. In this setting, our inferential procedure may 
be modified accordingly to appropriately take this into account. Stephens (1994) considers 
continuity-constrained switching linear regressions, and uses rejection sampling to draw from 
the conditional distribution of the changepoint 7 given the remaining model parameters. 
A more effective procedure, however, is to transform the data such that all but one of the 
regression coefficients may be integrated out: this variance reduction yields a Gibbs sampler 
analogous to Algorithm 1 which we detail below. In fact, it is possible to derive such an 
algorithm for any value of s > 0, though for ease of presentation we first describe the case 
s = 0, whereupon only continuity is enforced. 

This approach Jo variance reduction in the presence of srnoothness constraints proceeds 
as follows: define f3i = /3i and, without loss of generality, let /92 contain the last d2 elements 



1 D. Rudoy, S. G. Yuen, R. D. Howe and R J. Wolfe 

of /32, so that f3 = (/Jj^j/Jj) contains all p = di + d2 + 1 independent coefficients. Via tlie 
continuity constraint of (7), define a linear transformation oi (3 to (3 as follows: 

(I{dl + l)x{dl + l) 0(dl+l)xd2\ 
Clx(di+1) Oixrfs 1 > (11) 

0<i2X(dl + l) Id2Xd2 I 

where Imxm is the mxm identity matrix, Omxn is the mxn matrix of zeros, and Cix{di+i) — 
(1, a;^, a;^, . . . , x^^). The choice of which of the di+2 regression coefficients to select as the 
dependent variable is made without loss of generality, since a transformation similar to (11) 
can be defined for every such choice, as well as for the continuity constraint of (8). 

Since one of the regression coefficients is now a deterministic function of the others 
and the changcpoint, we place a prior directly on (3 rather than on (3. We assume /3i ^ 
J\f{fii,afA^^) and (32 ~ A/'(^2, 0-2-^2^)^ "^i^h A2 a (i2 x d2 diagonal positive definite matrix. 
Using the transformation Ty of (11), we obtain in analogy to (9) the full posterior 

• (|S^| |S|)-^ exp [-1 {{y - X,^)'S-i(y - X^P) + 0- /i)'S-iA(^ - ^i)}] , (12) 

where G j^nxp g^j^^^ analogy to the quantities (S,/i,A), we define S = 

diag(a2lrf,+i,aild,) G M^xp, fi = (K,/!^)' € W\ and A = diag(Ai, A2) G W^. 

The transformation T,,, makes it possible to integrate out the regression coefficients f3 
using standard manipulations. Indeed, introducing the terms = X!yH~'^Xj + S^^A G 
MP^P and = A'S~^jl + X'^'S~^y G MP^^ as before, we obtain the marginal posterior 

■ |S| exp {-i(y'S:7iy + jii'S-i/i - . (13) 

It is straightforward to generalize this notion to any s G {—1, 0, . . . , rfi + ^2}; a prior is put 
on di + (^2 — s + 1 regression coefficients, and a transformation T.^ analogous to (11) defined. 

As noted previously, the smoothness constraint of (7) introduces dependence among 
the pre- and post-changepoint regression coefficients. In contrast to the marginal posterior 
distribution of (9) derived earlier for the unconstrained case, enforcement of (7) precludes 
integrating out the associated noise variances and CTj. In the former case, the block- 
diagonal structure of X^ (and therefore of Aj) implies that the induced joint distribution 
on the variances is separable. However, in the latter case of (13), X^ is not block diagonal — 
owing to the action of — and hence neither is A^. Therefore, the variances af and cri 
are no longer conditionally independent, and their joint distribution does not take the form 
of known generalizations of the univariate Gamma distribution to the bivariate case (Yue 
ct al., 2001). Consequently, only the conditional distribution of the hyperparametcr bo is 
in standard form, and simulation from (13) proceeds with all other variables drawn using 
Metropolis-Hastings (MH) steps, as shown in Algorithm 2 below. In contrast to the case 
of Algorithm 1, where a mixture kernel was employed purely for reasons of computational 
efficiency, we emphasize here that such a move is in fact required to sample from the full 
support (l,n) of the changepoint index, otherwise mixing of the underlying chain is poor. 



Bayesian changepoint analysis for atomic force microscopy and soft materiai indentation 1 1 

As before, the mixture kernel consisted of a local random walk move and an independent 
global move drawing from a discrete distribution derived as a pointwise evaluation of (13) 
on the integers 1, 2, . . . , n. The coupling of noise variances suggests a joint MH move, but 
this requires specification of a proposal covariance; in simulations we found separate Normal 
random walk moves for ln(crf ) and ln(c72) to be adequate. 



Algorithm 2 Smoothness-constrained Gibbs sampler for changepoint estimation 

(a) Draw 7 ^ p{j\af,a2,bo, y; xp) according to (13) using a Metropolis- within-Gibbs step; 

(b) Draw ~ p(o'i I7, (t|, bo, y; xp) likewise; 

(c) Draw erf p(a-||7, af, bo, y, ip) likewise; 

(d) Draw 60 ~ p(&o|7, f^i, cri; 2/; V') = S{K,r]-^ +crf^ + fT^^). 



We re-emphasize that in our experience, integrating out the regression coefficients is 
essential in order to obtain a Gibbs sampler with favorable mixing properties. In particular, 
a sampler drawing for each parameter of (12) in the presence of smoothness constraints is 
severely restricted in its ability to explore the state space if one-coordinate-at-a-time updates 
are employed. When s continuous derivatives are enforced at the changepoint, the likelihood 
function becomes more nonlinear in changepoint Xj as s increases. Thus, as the induced 
constraint set becomes more nonlinear, local moves on the scale of the regression coefficients 
themselves will lie far from it — leading to small acceptance probabilities. To overcome these 
problems, one may design a high-dimensional MH move to update all the smoothness- 
constrained regression coefficients jointly with the changepoint; however, a unique move 
must be designed for each s to be considered. In contrast, integrating out the regression 
coefficients obviates this need by handling such constraints for all s simultaneously. 

5. Experimental validation 

In order to experimentally validate Algorithms 1 and 2 prior to their application in the 
setting of atomic force microscopy, we designed and performed two sets of special macro- 
scale indentation experiments in which precise contact point identification was made possible 
by the use of an impedance-measuring electrode mounted on the indenter. While this 
direct measurement procedure is precluded in the vast majority of biomaterials experiments 
involving cells and tissues, as such samples are siibmcrgcd in an aqueous solution, we were 
able to measure the true contact point for experiments involving respectively cantilever 
bending and silicone indentation, as described below. 

In addition, we conducted a number of simulation studies to characterize uncertainty in 
changepoint estimation as a function of various model parameters. On the basis of these 
simulation studies and subsequent experimental validation, we found both algorithms to 
be insensitive to the exact choice of hyperparamcters xp, and hence retained the following 
settings throughout: ^ = 0, A = 10~^ Jpxp, ao = 2, k = 1, and r] = 10~^, with p being the 
total number of regression coefficients. We set the pre-contact polynomial degree di = 1 
throughout, based on the discussion in Section 2 supporting the assumption of a linear pre- 
contact regime. When smoothness constraints were used, we retained identical parameter 
settings, and decremented the value of p appropriately. All posterior distributions were 
obtained by nmning the appropriate Gibbs samplers for 50 000 iterations, and discarding 
the first 5 000 samples. Convergence was assessed using standard methods confirming that 
increasing the number of Gibbs iterations did not appreciably change the resulting inference. 
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5. 1. Validation of changepoint inference via cantilever bending 

We first performed several trials of an experiment whereby a steel cantilever was bent 
through application of a downward uniaxial force. Here, the measured force F is expected 
to change linearly with displacement 5 = x — according to Hooko's Law, with 7 £ (1, n) 
representing the contact point index. Representative data shown in Figure 3 were obtained 
using a TestBench'^'^ Series system with a high-fidelity linear actuator (Bose Corporation 
EnduraTEC Systems Group, Minnetonka, Minnesota. USA), which moved a cylindrical 
indenter into a cantilevered piece of FSS-05/8-12 spring steel measuring approximately 
1.27 cm X 2.54 cm (carbon content 0.9-1.05%; Small Parts, Inc., Miramar, Florida, USA) 
at a speed of 10 mm/s. According to the impedance measurement technique described 
above, the contact point index 7 was determined to correspond to position index 48 of 
the indenter, corresponding to a contact point x-y G [—1.358, —1.262] mm. Because of the 
hardness of steel and the speed of the indenter, we did not make a smoothness assumption, 
and thus used the Gibbs sampler of Algorithm 1 to evaluate the efhcacy of our approach, 
with p = 4 based on the linear post-contact regime implied by Hooke's law. 

A full 100% of posterior values for 7 after a 10% burn- in portion took the value 48, 
indicating correct detection of the changepoint. Figure 3 shows the results of the curve 
fitting procedure,with similar results obtained for varying indenter speeds. The minimum 
mean-square error (MMSE) estimate of Young's modulus was determined to be 215.3 GPa, 
with an associated 95% posterior interval of [214.0,216.6]. By comparison, the range of 
values of Young's modulus for steel with similar carbon content is reported in the literature 
to be 210 ± 12.6 GPa (Kala and Kala, 2005). Despite its simple design, this experiment is 
not far removed from practice; a nearly identical procedure was employed by Wong et al. 
(1997) to probe the mechanical properties of silicone-carbide nanorods, with each nanorod 
pinned on a substrate and subjected to a bending force along the unpinned axis. 

5.2. Analysis of silicone indentation data 

While changes in slope at the contact point tend to be more pronounced in harder materials 
such as steel, the change in measured force is typically smoother in softer materials such 
as cells and tissues, making the contact point harder to detect. In earlier work (Yuen 
et al., 2007) we detailed an indentation experiment using a soft silicone sample (Aquaflex; 
Parker Laboratories, Fairfield, New Jersey, USA), chosen both for its similarity to human 
tissues used in material indentation studies (Chen et al., 1996), and because it enables direct 
contact point determination via an impedance-measuring electrode. 

Ten trials of this experiment were conducted, using a sample roughly 20 mm in depth, 
with a maximum indentation of approximately 8 mm; a typical force-displacement curve 
was shown earlier in Figure 2. A hemispherical metal indenter of radius R = 87.5 mm com- 
pressed the sample at 10 mm/s, and the resulting forces were measured approximately every 
10 /im to yield ~ 960 force-displacement data points. Both the imconstraincd (s = — 1) and 
continuity-constrained (s = 0) models were fitted in this setting, by way of Algorithms 1 
and 2 respectively. For a spherically-tipped indenter, the Hertzian model of (2) indicates 
that force is proportional to (x — x-y)^/^, and hence we employed a post-contact design ma- 
trix with only the fractional power 3/2. In this regime, we have that p = 4 for Algorithm 1 
and p = 3 for Algorithm 2. Since the initial distance from the indenter to the sample was 
approximately known, a imiform prior on 7 G [125, 250] was assumed. 

For each of the ten data sets collected, the first n = 450 data points were taken to 
represent a conservative estimate of operation within the small-deformation regime, and 
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Fig. 3. Inference for the cantilever experiment of Section 5.1. The force-displacement data and 
posterior mean reconstructions obtained via Algorithm 1 are shown (top left) with a close-up view 
near the changepoint (bottom left) and the associated residuals (top right). The marginal posterior of 
Young's modulus (bottom right) is shown with its mean (black line) and 95% interval (grey lines). 



the Gibbs samplers of Algorithms 1 and 2 were each run on these data. Results for both 
the cases are summarized in Table 1, along with the experimentally-determined contact 
point, which varied from trial to trial due to viscoclastic effects. Indeed, given that the 
spacing between data points is 0.01 mm on average, it may be deduced from Table 1 that 
the average error of 0.8-1% across trials corresponds to 8-10 sampled data points. 

Marginal contact point posterior distributions for both the unconstrained and continuity- 
constrained cases are summarized, pairwise by experiment, in the box plots of Figure 4. 
These are seen to be notably more diffuse in the former case (left) than the latter (right); 
this is consistent with the softness of the silicone sample under study, which makes it 
difficult to reject a priori the possibility of continuity at the changepoint. Absent this 
assumption, the experimentally determined contact point lies within the 95% posterior 
interval for each of the ten trials. Once a continuity constraint is imposed, the marginal 
posterior distributions of the changepoints tighten noticeably; this is consistent with our 
expectation that constraining the model reduces the number of high-likelihood fits. 

A more subtle point can also be deduced from the slight yet consistent left-shift across 
95% posterior intervals under the continuity-constrained regime relative to the uncon- 
strained case. Observe the fourth row of Table 1, which reports the values obtained upon 
extrapolating pre- and post-contact MMSE curve fits to their meeting point in this latter 
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Table 1. Contact point estimates and associated errors (mm, %) for the ten silicone trials of Sec- 
tion 5.2, based on unconstrained (U) and continuity-constrained (C) models, with force-response data 
sampled every 0.01 mm on average. Young's modulus estimates E and 95% posterior intervals for 
the unconstrained case are also shown, along with averages over all ten trials where appropriate. 



Trial No. 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


Avg. 


Truth (a;^) 


5.52 


5.49 


5.47 


5.44 


5.42 


5.38 


5.39 


5.48 


5.42 


5.38 




MMSE (U) 


5.42 


5.44 


5.42 


5.43 


5.47 


5.39 


5.46 


5.35 


5.41 


5.39 




MMSE (C) 


5.40 


5.43 


5.41 


5.39 


5.46 


5.37 


5.41 


5.34 


5.34 


5.39 




Extrap. (U) 


5.39 


5.42 


5.40 


5.37 


5.42 


5.36 


5.40 


5.32 


5.34 


5.38 




% Err. (U) 


-1.78 
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Fig. 4. Box plots of contact point marginal posterior distributions in each of ten silicone indentation 
trials, shown side-by-side for the unconstrained (left) and continuity-constrained (right) models, with 
centers of the grey squares indicating true contact point values. Note the consistent decrease in 
posterior variance and slight downward shift of the posterior under the assumption of continuity. 



case. The increase in force after contact implies that these values will always lie below the 
directly inferred contact point, as indeed they do. Nevertheless, enforcing continuity results 
in only a small increase in overall contact point estimation error, and produces what practi- 
tioners might judge to be more physically feasible curve fits. While independent verification 
of the Young's modulus is not available for the silicone sample used in our study, the quality 
of our contact point estimates relative to a known ground truth leads to high confidence in 
the inferred values of Young's modulus. 

Overall, the ability to obtain inferential results and accompanying uncertainty quan- 
tification, as exemplified by the cantilever bending and silicone indentation experiments. 
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represents a significant improvement upon current methods; a more in-depth comparison 
to the method of Costa et al. (2006) is provided in our earlier work (Yuen et al., 2007). 
In particular, the latter experiment demonstrates that reliable estimates of soft material 
properties can be obtained even in the presence of measurement error — a regime applicable 
to many AFM studies of cells and tissues, as we now describe. 

6. Inference in the setting of atomic force microscopy 

Having validated our algorithms on a macroscopic scale, we turn to analyzing cellular bio- 
materials data collected using atomic force microscopy techniques. In contrast to our earlier 
examples, no direct experimental verification is available in this case, though we note that 
our resultant estimates of Young's modulus are considered plausible by experimentalists 
(Socrate and Suresh Labs, personal communications, 2008). 

It is widely believed that cell biomechanics can shed light on various diseases of import, 
and hence a key research objective following the advent of AFM technology has been to 
analyze quantitatively the mechanical properties of various cell types. Indeed, numerous 
papers over the past decade have linked stiffness and other related mechanical properties 
to cell malfunction and death (Costa, 2004); for example, the ability of cardiac myocytes 
in heart muscle tissue to contract is intimately linked to their cytoskeletal structure and 
its influence on cellular mechanical response (Lieber et al., 2004). Here, we are similarly 
motivated to understand the stiffness properties of neuronal and red blood cells — currently 
a topic of intensive research in the biomechanics and bioengineering communities. 

6. 1. Indentation study of embryonic mouse cortical neurons 

We first analyzed ex- vivo live mouse neurons, submerged in cell culture medium and repeat- 
edly indented by an AFM (Asylum Research, Santa Barbara, California, USA) equipped 
with a spherically-tipped probe. The indentation of each neuron (Socrate Lab, Mas- 
sachusetts Institute of Technology) yielded approximately 700 force measurements, of which 
the first 500 were used in subsequent analysis in order to stay within the small- deformation 
regime. Such data arc of wide interest to ncuroscientists and engineers, as traumatic damage 
to neurons is hypothesized to be related to their mechanical properties. 

As a spherical probe tip was used for indenting each cell, we employed a post-contact 
design matrix with only the fractional power 3/2, as in the case of our earlier silicone ex- 
ample, and the continuity-constrained sampler of Algorithm 2. The results of a typical 
trial are shown in Figure 5; the pre- and post-contact residuals were observed to be white. 
The primary inferential quantity of interest in such cases is Young's modulus, the princi- 
pal characterization of cell material stiffness introduced in Section 2.2. According to the 
Hertzian model of (2) for a spherical indenter, the regression coefficient corresponding to the 
(a; — a;^)^/^ term of the fitted model is proportional to Young's modulus, with the constant 
of proportionality a function of the given radius i? = 10 /um and Poisson's ratio u = 0.5. 
Thus, we can obtain the posterior distribution of Young's modulus by appropriately scaling 
the distribution of this regression coefficient, as shown in the bottom left panel of Figure 5. 
We report the MMSE estimate of Young's modulus as E = 530.4 Pa, and the correspond- 
ing 95% posterior interval as [518.8,544.1] Pa. This estimate is in reasonable agreement 
with those previously reported for similar neurons (Lu et al., 2006; Elkin et al., 2007), 
with variability due primarily to differences in indentation speeds (Socrate Lab, personal 
communication, 2008). 
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Fig. 5. Data collected during an AFM indentation of a mouse neuron together with the posterior mean 

estimate of the underlying regressions (top left) and the induced residuals (top right). Also shown 
are marginal posterior distributions of Young's modulus (bottom left) and contact point (bottom right), 
each overlaid with the posterior mean (black line) and 95% posterior interval (dashed grey lines). 



6.2. Indentation study of red blood cells 

The mechanics of red blood cells have also been extensively studied using atomic force 
microscopy (Radmachcr ct al.. 1996). In this vein, we next analyzed data from ex-vivo live 
human erythrocytes (red blood cells, Suresh Lab, Massachusetts Institute of Technology) 
submerged in a cell culture mcdimn, indented by an AFM (Asylum Research) equipped with 
a pyramidally-tipped probe. The final 800 of approximately 8 500 data points were discarded 
prior to analysis, as they clearly lay outside the small- deformation regime. Relative to the 
neuron AFM data considered in Section 6.1, the sampling rate of resistive force here is high, 
and consequently it is feasible to enforce continuity (s = 0) at the changepoint. Algorithm 2 
was therefore employed, with results from a typical trial shown in Figure 6. 

In the case at hand, the regression coefficient /3i2 corresponding to (x — x-yY is propor- 
tional to the Young's modulus as per (2), with the constant of proportionality depending 
on the indenter tip angle 20 = 70° and Poisson's ratio v = 0.5, and, as we discuss below, 
no linear post-contact regression term was included. The inferred distribution of Young's 
modulus may be obtained by appropriately transforming the marginal posterior of /3i2, as 
detailed in Section 3.1, and is shown in the bottom left panel of Figure 6. The resultant 
MMSE estimate of E = 25.3 Pa and corresponding posterior interval of [16.0, 34.9] Pa were 
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Fig. 6. Data collected during an AFM indentation of a human red blood cell together with the posterior 
mean estimate of the underlying regressions (top left) and the induced residuals (top right). Also 
shown are marginal posterior distributions Young's modulus (bottom left) and of al and al (bottom 
right), the former overlaid with the posterior mean and 95% posterior interval. 



confirmed to be eonsistent with various experimental assumptions (Suresh Lab, personal 
communication, 2008). Further, the pre- and post-contact error variances are determined 
to be unequal, as shown in the bottom right panel of Figure 6. 

Note that in the absence of a post-contact drift, enforcing continuous differentiability 
at the changepoint (s = 1) constrains the pre-contact linear fit to have zero slope, and is 
inconsistent with the pre-contact drift clearly visible in the top left panel of Figure 6. On the 
other hand, if the post-contact polynomial were to include a linear drift term, then enforcing 
continuous differentiability would imply that the pre- and post-contact drifts are identical. 
Though this is appealing from a modeling viewpoint, as it eliminates an additional free 
parameter, practitioners lack evidence for such an equality. Moreover, in our experiments 
its inclusion had no appreciable effect on the inference of Young's modulus, and so we did 
not include a post-contact drift term in our final analysis of these data. 

As in the case of our earlier experiments, we compared our approach to the likelihood 
method of Costa et al. (2006), which yielded an estimate of shifted to the right by more 
than 1000 data points relative to that shown in Figure 6. The corresponding estimate 
of Yoimg's modulus in turn was found to be 34.8 Pa — an increase of 37.5% relative to the 
MMSE point estimate of 25.3 Pa, and close to the upper boundary of our estimated posterior 
interval. While in this experimental setting one cannot conclude that either estimate is 
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superior to the other, wc note that the difference between them can be attributed in part to 
our model's incorporation of differing pre- and post-contact error variances, and smoothness 
constraints. Thus, one may view our inferential procedures as both a formalization and an 
extension of earlier likelihood-based approaches developed by practitioners, enabling both 
robust, automated fitting procedures and explicit uncertainty quantification. 

7. Discussion 

In this article we have posed the first rigorous formulation of — and solution to — the key 
inferential problems arising in a wide variety of material indentation systems and studies. 
In particular, practitioners in the materials science community to date have lacked accurate, 
robust, and automated tools for the estimation of mechanical properties of soft materials 
at either macro- or micro-scales (Lin et al., 2007a; Crick and Yin, 2007). A principal 
strength of our approach is its applicability to the analysis of biomaterials data obtained 
by indenting cells and tissues using atomic force microscopy; contact point determination 
is even more difficult in this setting, owing to the gradual change of measured resistive 
force that is a hallmark of soft materials. The Bayesian switching polynomial regression 
model and associated inferential procedures we have proposed provide a means both to 
determine the point at which the indenting probe comes into contact with the sample, and 
to estimate the corresponding material properties such as Young's modulus. In turn, its 
careful characterization holds open the eventual possibility of new biomechanical testing 
procedures for disease (Costa, 2004). 

Our parametric approach is strongly motivated by — and exploits to full advantage — 
the Hertzian models governing the physical behavior of linear-elastic materials undergoing 
small deformations. The Bayesian paradigm not only enables uncertainty quantification, 
crucial in applications, but also allows for the natural incorporation of physically- motivated 
smoothness constraints at the changepoint. Inference is realized through application of 
carefully designed Markov chain Monte Carlo methods together with classical variance re- 
duction techniques. The resultant algorithms have been shown here to be both statistically 
and computationally efficient as well as robust to choice of hyperparameters over a wide 
range of examples, and are available online for use by practitioners. Indeed, the direct 
applicability of our methods precludes any need for data pre-processing prior to analysis. 

Outside of the linear-elastic materials we consider here, it is of interest to apply the 
methodology to viscoelastic materials (e.g., biopolymers) which return to their pre-contact 
state slowly over time (Lin et al., 2007b). In such cases, the amount of induced deformation 
depends not only on the indenter geometry, but also on the rate of indentation. Another 
extension is to incorporate multiple spatially distributed changepoints into our model, a 
key construct when atomic force microscopes are used to indent repeatedly a sample in 
order to characterize cell stiffness as a function of surface location (Geisse, 2009). Finally, a 
sequential estimation scheme could be of great use in surgical robotics applications, where 
contact point determination plays a key role in enabling tactile sensing — a subject of current 
study by the authors. 
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